rm(list=ls())
setwd("C:/Users/bschneer/Documents/GOV1/GOV 2001/Replication/CEM")
library(cem)
library(mvtnorm)
library(modeest)
library(Matching)
sims<-1000000
set.seed(02138)

files<-c("ada2","capital2","federalism2","sexdisc2","sexhar2","takings2","titlevii2")
diagnostic<-c()
bal<-c()
holder<-as.data.frame(matrix(0,nrow=length(files),ncol=7))
names(holder)<-c("Data","ATE","Lower","Upper","L1 (Pre-CEM)","L1 (Post-CEM)","T in Sample After Match")

beta<-vector("list",length(holder))
var<-vector("list",length(holder))
Z<-vector("list",length(holder))
W<-vector("list",length(holder))

load("genderjudging.RData")

for (i in 1:length(files)){
match<-get(files[i])[,c("judge.vote","treatment","circuit","dec.year","jcs","year.birth","minority.judge","lower.dir","abs.jcs","jud.experience")]
match["circuit"]<-as.factor(match[,"circuit"])
match<-match[match$dec.year<2003,]

imb1<-imbalance(match$treatment,match,drop=c("judge.vote","treatment"))

dec.year_cut<-list(c(1998, 1999, 2000, 2001, 2002),
c(1995, 1996.75, 1998.5, 2000.25, 2002),
c(1995, 1996.17, 1997.33, 1998.5, 1999.67, 2000.83, 2002),
c(1995, 1996.17, 1997.33, 1998.5, 1999.67, 2000.83, 2002),
c(1995, 1996.33, 1997.67, 1999, 2000.33, 2001.67, 2003),
c(1979, 1983.6, 1988.2, 1992.8, 1997.4, 2002),
c(1995, 1996.75, 1998.5, 2000.25, 2002))

jcs_cut<-list(c(-0.78, -0.49, -0.2, 0.09, 0.38, 0.66),
c(-0.62, -0.32, -0.02, 0.28, 0.58),
c(-0.62, -0.31, 0.01, 0.32, 0.64),
c(-0.79, -0.33, 0.12, 0.58),
c(-0.73, -0.4, -0.07, 0.25, 0.58),
c(-0.66, -0.41, -0.15, 0.1, 0.36, 0.61),
c(-0.78, -0.54, -0.31, -0.07, 0.17, 0.4, 0.64))

year.birth_cut<-list(c(1905, 1916.2, 1927.4, 1938.6, 1949.8, 1961),
c(1910, 1924.67, 1939.33, 1954),
c(1901, 1918.67, 1936.33, 1954),
c(1905, 1930, 1955),
c(1906, 1922.67, 1939.33, 1956),
c(1897, 1906.5, 1916, 1925.5, 1935, 1944.5, 1954),
c(1906, 1964))

minority.judge_cut<-list(c(0, .5),
c(0, .5),
c(0, .5),
c(0, 1),
c(0, 1),
c(0, 0.5),
c(0, .5))

lower.dir_cut<-list(c(0, 0.5),
c(0, 0.5),
c(0, 0.5),
c(0, 0.5),
c(0,.5),
c(0, 0.5),
c(0,.5))

abs.jcs_cut<-list(c(0, 0.18, 0.36, 0.53, 0.71, 0.89, 1.07),
c(0, 0.3, 0.6, 0.9),
c(0, 1.07),
c(0, 0.53, 1.07),
c(0, 0.2),
c(0, 0.56, 1.12),
c(0, 0.2))

jud.experience_cut<-list(c(0, .5),
c(0, 1),
c(0, 1),
c(0, 1),
c(0, 1),
c(0,1),
c(0,1))

circuit_gps<-list("1","2","3","4","5","6","7","8","9","10","11","DC")

cp<-list(dec.year=dec.year_cut[[i]],jcs=jcs_cut[[i]],year.birth=year.birth_cut[[i]],minority.judge=minority.judge_cut[[i]],lower.dir=lower.dir_cut[[i]],abs.jcs=abs.jcs_cut[[i]],jud.experience=jud.experience_cut[[i]])
gps<-list(circuit=circuit_gps)

cem.match<-cem(treatment="treatment",data=match,cutpoints=cp,grouping=gps,drop=c("judge.vote"))
cem.match
imb2<-imbalance(match$treatment,match,drop=c("judge.vote","treatment"),weights=cem.match$w)

m1<-MatchBalance(treatment~jcs + year.birth + minority.judge + lower.dir + abs.jcs + dec.year, data=match)
m2<-MatchBalance(treatment~jcs + year.birth + minority.judge + lower.dir + abs.jcs + dec.year, data=match, weights=cem.match$w)

before1<-c()
after1<-c()
for (j in 1:6) before1<-cbind(before1,m1$BeforeMatching[[j]][3:6])
for (j in 1:6) after1<-cbind(after1,m2$BeforeMatching[[j]][3:6])

before1<-data.frame("before",before1)
after1<-data.frame("after",after1)

names(before1)<-c("t","jcs","year.birth","minority.judge","lower.dir","abs.jcs","dec.year")
names(after1)<-c("t","jcs","year.birth","minority.judge","lower.dir","abs.jcs","dec.year")
mns<-rbind(before1,after1)

bal[[i]]<-mns

#diagnostic[[i]]<-imbspace(cem.match,match,M=250,depth=3,plot=F)

match$dec.era<-as.numeric(match$dec.year<1988)
match$dec.era[match$dec.year>=1988 & match$dec.year <1992]<-2
match$dec.era[match$dec.year>=1992 & match$dec.year <2000]<-3
match$dec.era[match$dec.year>=2000]<-4

spec1<-formula(judge.vote ~ treatment)
spec2<-formula(judge.vote ~ treatment + jcs + year.birth + minority.judge + lower.dir + abs.jcs + as.factor(dec.year))
spec2b<-formula(judge.vote ~ treatment + jcs + year.birth + minority.judge + lower.dir + abs.jcs + as.factor(dec.era))


model<-glm(spec2,binomial(logit),data=match,weights=cem.match$w)
if(i==6){
model<-glm(spec2b,binomial(logit),data=match,weights=cem.match$w)
}
summary(model)

mean(fitted(model))*mean(1-fitted(model))*coefficients(model)[2]

beta[[i]]<-model$coefficients
var[[i]]<-summary(model)$cov.unscaled
Z[[i]]<-as.data.frame(model.matrix(model))
W[[i]]<-cem.match$w

invlogit <- function(x){
  1 / (1 + exp(-x))
}
ate1 <- function(model, sims){
  data        <- model$data
  beta        <- coef(model)
  V           <- summary(model)$cov.scaled
  beta.draw   <- rmvnorm(sims, mean=beta, sigma=V)
  pr.men      <- matrix(0, nrow(as.matrix(data)), sims)
  pr.women    <- matrix(0, nrow(data), sims)
  X.men       <- X.women <- model.matrix(model)
  X.men[,2]   <- 0
  X.women[,2] <- 1

  diff.holder <- matrix(0, nrow(data), sims)

  for(i in 1:nrow(data)){
    p0   <- invlogit(X.men[i,] %*% t(beta.draw))
    p1   <- invlogit(X.women[i,] %*% t(beta.draw))
    diff <- p1 - p0
    
    diff.holder[i,] <- diff
  }
  return(apply(diff.holder, 2, mean))
}

############################
##Ave Treatment of Treated##
############################

sim<-ate1(model,1000)
holder[i,]<-cbind(files[i],mean(sim),quantile(sim,.025),quantile(sim,.975),imb1$L1$L1,imb2$L1$L1,cem.match$tab[2,2]/cem.match$tab[1,2])
}
holder

save(bal,file="balance_cem.RData")
save(holder,file="estimates.RData")

load("pscore_L1.RData")
comparison<-cbind(holder,L1)
names(comparison)[8]<-"Propensity Score L1"

write.csv(comparison,file="pscore_compare.csv")

#pdf(paste("diagnostic",i,".pdf",sep=""),width=11,height=8.5)
#par(mfrow=c(2,3),oma=c(0,0,0,0))
#for (i in 1:6) plot(diagnostic[[i]],data=match,explore=F)
#dev.off()

